home *** CD-ROM | disk | FTP | other *** search
/ World of Education / World of Education.iso / world_e / esptest.zip / ESPPROB.PAS < prev    next >
Pascal/Delphi Source File  |  1992-03-05  |  12KB  |  347 lines
  1. UNIT ESPPROB;
  2. {**********************************************************}
  3. {**                                                      **}
  4. {**      A Unit to compute probability  for ESPTEST      **}
  5. {**      Copyright 1991 Phil Mosier                      **}
  6. {**                                                      **}
  7. {**      Used with permision                             **}
  8. {**      From BIO-STATISTICAL MICROCOMPUTING IN PASCAL   **}
  9. {**      BY Klotz and Meyer                              **}
  10. {**      (c) 1985 Rowman & Allanheld Publishers          **}
  11. {**                                                      **}
  12. {**********************************************************}
  13. Interface
  14. Uses Graph, Crt;
  15. Procedure BINOMIAL(Var PROB_B: Real;
  16.                        TRIALS: Word;
  17.                        HITS:   Word);
  18.  
  19. Procedure CHISQUARE(Var PROB_C: Real;
  20.                         CHI_SQ: Real;
  21.                         DEGREE_FREEDOM: Real);
  22.  
  23. Implementation
  24. Procedure BINOMIAL(Var PROB_B: Real;
  25.                        TRIALS: Word;
  26.                        HITS:   Word);
  27.  
  28. {**********************************************************}
  29. {**    From BIO-STATISTICAL MICROCOMPUTING IN PASCAL     **}
  30. {**    BY Klotz and Meyer                                **}
  31. {**********************************************************}
  32.  
  33. Var
  34. X,XL,A,B,N,P,PL : Real;
  35. {*********************************************************}
  36. {**                                                     **}
  37. {**                                                     **}
  38. {*********************************************************}
  39. Function BINET(A:Real): Real; { Benet's Function }
  40.    Var
  41.    S: Real;
  42.    Begin {BINET}
  43.    {*************************  BINET  ********************}
  44.       If A < 1.5 Then
  45.          BINET := 8.1061467E-2
  46.         {*************************************************}
  47.         {**      Case Binet(1) = 1.5*ln(2pi)            **}
  48.         {*************************************************}
  49.      Else
  50.         {*************************************************}
  51.         {**    Continued Fraction                       **}
  52.         {*************************************************}
  53.         Begin
  54.            S := A + 1.5174736 /( A + 2.269489);
  55.            S := A + ( 195 /371 ) /( A + 22999 /22737 /S);
  56.            BINET := 1 /12 /( A + ( 1 /30 ) /( A + ( 52 /210 ) /S));
  57.         End;
  58. End; {BINET}
  59.  
  60. Function G(U : Real) : Real;
  61. {*********************************************************}
  62. {**   LN(1+U)-U With Accuracy                           **}
  63. {*********************************************************}
  64. Var
  65. R: Real;
  66. Begin {G}
  67.    If ((U < -0.5) or ( U > 1.0)) Then
  68.       G := Ln(1 + U) - U
  69.    Else
  70.    {*************************************************}
  71.    {**   continued fraction                        **}
  72.    {*************************************************}
  73.    Begin
  74.       R:= 6 + 16 * U /(7 + 9 * U /(8 + 25 * U /(9 + 16 * U /10)));
  75.       G := -U * U /(2 + 4 * U /(3 + U /(4 + 9 * U /(5 +4 * U /R ))));
  76.    End;
  77. End;  { G}
  78.  
  79. Function ICBETA(X,XL,A,B : Real) : Real;
  80. {*********************************************************}
  81. {**  Program for the incomplete beta distribution.      **}
  82. {**  The algorithm of H.E. Soper (1921), Tracts for     **}
  83. {**  Computer No. VII, Karl Person Ed., Cambridge Univ. **}
  84. {** Press, is used with a,b incremented if either < 1.  **}
  85. {*********************************************************}
  86. Var
  87. C,FXAB : Real;
  88.    Function F(X,XL,A,B : Real) : Real;
  89.    Var
  90.    S,U,V,GU,GV : Real;
  91.       Begin {F}
  92.       U := X * (A + B) /A - 1;
  93.       V := XL * A /B - X;
  94.       If U < -0.9 Then GU := Ln(X * (A + B) /A) - U
  95.       Else GU := G(U);
  96.       If V < -0.9 Then GV := Ln(XL * (A + B) /B) - V
  97.       Else GV := G(V);
  98.       S := 0.5 * Ln(A * B /(( A + B) * 6.2831853));
  99.       S := S + A * GU + B * GV + BINET(A + B) - BINET(A) - BINET(B) - Ln(A);
  100.       If S < - 85.19
  101.          {************************************************}
  102.          {**  Ln(IE-37)                                 **}
  103.          {************************************************}
  104.          Then F := 0
  105.          Else F := Exp(S);
  106.    End; {F}
  107.   
  108.    Function SOPER(X,XL,A,B : Real) : Real;
  109.    {********************************************************}
  110.    {**                                                    **}
  111.    {********************************************************}
  112.    Const
  113.    EPS = 1.0E-7;
  114.    Var
  115.    V,AI,CK,SUM,ABL,AS,A12,B2,AB4,R : Real;
  116.    S,SMAX,I,K : Integer;
  117.    Begin { SOPER  }
  118.    {*************************  soper  ********************}
  119.       SUM := 0;
  120.       I := 0;
  121.       AI := 1;
  122.       V := X /XL;
  123.       R := B + XL * (A + B);
  124.       SMAX := MAXINT - 1;
  125.       If R < SMAX Then
  126.          S := TRUNC(R)
  127.       Else S := SMAX;
  128.       If S > 0 Then
  129.          Repeat
  130.          {****************************************************}
  131.          {**    Raise a Decreas b                           **}
  132.          {****************************************************}
  133.             SUM := SUM + AI;
  134.             I := I + 1;
  135.             AI := AI * V * (B - I ) /(A + I);
  136.          Until ( Abs(AI /SUM) < EPS) or ( I >= S);
  137.       If ( I = S ) Then
  138.          Begin
  139.          {**************************************************}
  140.          {**    Raise a                                   **}
  141.          {**************************************************}
  142.              K := 0;
  143.              CK := XL * AI;
  144.              ABL := A + B - 1;
  145.              AS := A + S;
  146.              Repeat
  147.                 SUM := SUM + CK;
  148.                 K := K + 1;
  149.                 CK := X * CK * (ABL + K) /(AS + K);
  150.             Until (ABS(CK /SUM) < EPS);
  151.       End; { raise a}
  152.       A12 := (A + 1) * (A + 2);
  153.       B2 := B * (B + 1);
  154.       AB4 := (A + B) * (A + B + 1) * (A + B + 2) * (A + B + 3);
  155.       SOPER := SUM * F(X, XL, A + 2, B + 2) * A12 *
  156.           B2 /AB4 /SQR(X * XL) /XL;
  157.    End; {soper}
  158.  
  159. Begin
  160. {****************************  ICBETA  *******************}
  161. If X <= 0 Then
  162.    ICBETA := 0
  163. Else
  164.    If X >= 1.0 Then
  165.       ICBETA := 1.0
  166.    Else
  167.       If (A > 1) And (B > 1) Then
  168.          If X <= (A /(A - B)) Then
  169.             ICBETA := SOPER(X, XL, A, B)
  170.          Else ICBETA := 1.0 - SOPER(XL, X, B, A)
  171.       Else
  172.          Begin
  173.          {********************************************}
  174.          {**    Increment a,b                       **}
  175.          {********************************************}
  176.          C:= A * (A + 1) /(A + B) * (A + 2) /( A + B + 1) *
  177.             B /(A + B + 2) * (B + 1) /(A + B + 3);
  178.          FXAB := F(X, XL, A + 2, B + 2) * C /Sqr(X * XL);
  179.          ICBETA := FXAB * (1 - X * (A + B) /B) /A +
  180.             ICBETA(X, XL, A + 1, B + 1);
  181.          End;
  182.       {**********************************************}
  183.       {**     Increment a,b                        **}
  184.       {**********************************************}
  185. End;  { ICBETA}
  186.  
  187. Begin
  188. {***************************  Probability  ***************}
  189.       N := TRIALS;
  190.       P  := 0.2;
  191.       PL := 0.8;
  192.       X := HITS;
  193.       PROB_B := ICBETA(PL, P, N - X, X + 1);
  194. End;          {probabilty}
  195.  
  196. Procedure CHISQUARE(Var PROB_C : Real;
  197.                         CHI_SQ : Real;
  198.                         DEGREE_FREEDOM : Real);
  199. Var
  200. X,N : real;
  201. A : Real;
  202. TEMP : String[18];
  203. {*********************************************************}
  204. {**                                                     **}
  205. {**     Procedure CHISQUARE                             **}
  206. {**                                                     **}
  207. {*********************************************************}
  208. Function BINET(A:real) :Real; { Benet's Function }
  209.    Var
  210.    S : Real;
  211.    Begin {BINET}
  212.    {*************************  binet  ********************}
  213.       If A < 1.5 Then
  214.          BINET := 8.1061467E-2
  215.         {*************************************************}
  216.         {**      Case BINET(1) = 1.5*ln(2pi)            **}
  217.         {*************************************************}
  218.      Else
  219.         {*************************************************}
  220.         {**    Continued Fraction                       **}
  221.         {*************************************************}
  222.         Begin
  223.            S := A + 1.5174736 /(A + 2.269489);
  224.            S := A + (195 /371) /(A + 22999 /22737 /S);
  225.            BINET := 1 /12 /(A + (1 /30) /(A + (52 /210) /S));
  226.         End;
  227. End; {BINET}
  228.  
  229. Function G(U : Real) : Real;
  230. {*********************************************************}
  231. {**   ln(1+u)-u with accuracy                           **}
  232. {*********************************************************}
  233. Var         R  : Real;
  234. Begin {G}
  235.    {  book corrected)
  236.    If (U < -0.5) Or (U > 1.0) Then
  237.       G := Ln(1 + U) - U
  238.    Else
  239.    {*************************************************}
  240.    {**   continued fraction                        **}
  241.    {*************************************************}
  242.    Begin
  243.       R := 6 + 16 * U /(7 + 9 * U /(8 + 25 * U /(9 + 16 * U /10)));
  244.       G := -U * U /(2 + 4 * U /(3 + U /(4 + 9 * U /(5 + 4 * U /R))));
  245.    End;
  246. End;  { G}
  247.  
  248. {************************************************}
  249. {**                                            **}
  250. {**    Function LNGAM2                         **}
  251. {**                                            **}
  252. {************************************************}
  253. Function LNGAM2(A:Real) : Real ; { LN(GAMMA(2 + A)) }
  254. Var
  255. K : Integer;
  256. S : Real ;
  257. C : Array[1..9] Of Real; { C[1]=1-GAMMA, C[K]=(-1**K)*(ZETA(k)-1)/k}
  258. {******************  LNGAM2  *******************}
  259. Begin {LNGAM2}
  260.    C[1] :=  0.4227843; C[2] :=  0.3224670; C[3] := -0.0673523;
  261.    C[4] := -0.0205808; C[5] := -0.0073856; C[6] :=  0.0028905;
  262.    C[7] := -0.0011928; C[8] :=  0.0005097; C[9] := -0.0002232;
  263.    S := C[9] * A;
  264.    For K := 8 DownTo 1 Do S:= C[K] + A * S;
  265.       LNGAM2 := A * S;
  266. End; {LNGAM2}
  267.  
  268. {************************************************}
  269. {**                                            **}
  270. {**    Function ICGAMMA                        **}
  271. {**                                            **}
  272. {************************************************}
  273. Function ICGAMMA(A,X : Real) : Real; {Incomplete GAMMA C.D.F.}
  274. Var
  275. K: Integer;
  276. D,S,C,UK,AK : Real;
  277. {************************************************}
  278. {**                                            **}
  279. {**    Function LNXF                           **}
  280. {**                                            **}
  281. {************************************************}
  282.    Function LNXF(X,A : Real) : Real; {Ln(X*F(X,A))}
  283.    Var
  284.    U,GU : Real;
  285.    {********************  LNXF   ****************}
  286.    Begin {LNXF}
  287.       U := (X - A) /A;
  288.       If U < -0.9 Then GU := Ln(X /A) - U
  289.       Else GU := G(U);
  290.       If A < 0.5 Then LNXF := A * GU + G(A) + (A + 1) * LN(A) - LNGAM2(A)
  291.       Else
  292.          If A<1.5 Then
  293.             LNXF := A * GU + G(A - 1) + (A * LN(A) - 1) - LNGAM2(A - 1)
  294.          Else {case a >= 1.5}
  295.             LNXF := A * GU - BINET(A) - 0.5 * Ln(6.2831853 /A);
  296.    End; {LNXF}
  297.  
  298. {********************  ICGAMMA  **************}
  299. Begin {ICGAMMA}
  300.    If X <= 0 Then ICGAMMA := 0
  301.    Else
  302.       Begin {nonzero case}
  303.          D := Exp(LNXF(X,A));
  304.          If (X <= A) Or (X <= 1) Then
  305.             Begin {series}
  306.                AK := A;
  307.                UK := 1;
  308.                S := 1;
  309.                While (Abs(UK /S) > 1E-7) Do
  310.                   Begin
  311.                      AK := AK +1;
  312.                      UK :=UK * X /AK;
  313.                      S := S + UK;
  314.                   End;
  315.                ICGAMMA := S * D /A;
  316.             End {series}
  317.          Else
  318.             Begin {continued fraction}
  319.                C:= 0;
  320.                For K := 30 DownTo 1 Do
  321.                   Begin
  322.                      C := K /(X + C);
  323.                      C := (K - A) /(1 + C);
  324.                   End;
  325.                C:= 1 /(X + C);
  326.                ICGAMMA := 1 - C * D;
  327.             End;  {continued fraction}
  328.       End;  {non zero case}
  329. End; {ICGAMMA}
  330.  
  331. Begin { CHISQUARE}
  332. {**************************  CHISQUARE  *****************}
  333. N :=  DEGREE_FREEDOM;
  334. X :=  CHI_SQ;
  335. If N < 0.5 Then Halt(0)
  336. Else
  337.    Begin
  338.       If X < 0 Then PROB_C := 0
  339.       Else PROB_C := ICGAMMA(N /2.0,X /2.0);
  340.    End;
  341. End; { CHISQUARE }
  342.  
  343. End.    { Turbo Pascal Unit  --  ESP_PROB}
  344.  
  345.                                     
  346.  
  347.